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Abstract 

A new Monte Carlo scheme based on the system of Tsahis's gener- 
ahzed statistical mechanics is applied to a simple double well potential 
to calculate the canonical thermal average of potential energy. Al- 
though we observed serious quasi-ergodicity when using the standard 
Metropolis Monte Carlo algorithm, this problem is largely reduced by 
the use of the new Monte Carlo algorithm. Therefore the ergodicity is 
guaranteed even for short Monte Carlo steps if we use this new canon- 
ical Monte Carlo scheme. 
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1 Introduction 



The ergodic hypothesis is fundamental to statistical mechanics. This hy- 
pothesis states that the time average of an observable event equals the 
phase-space average. In practical application, however, problems can arise 
in various types of simulations if the system must overcome high energy 
barriers to reach other regions of phase space. In that case, the length 
of a simulation needed in order to obtain enough statistical samples of all 
regions of phase space may be extremely long. In the Monte Carlo simu- 
lation, the errors arise as a consequence of the finite length of the Monte 
Carlo walk. This error can be serious in canonical Monte Carlo sampling, 
especially at low temperatures. This problem, referred as "quasi-ergodicity" 
by Valleau and Whittington |jl|, appears even in the simplest double well 
problem, where the two wells are separated by a large barrier . 

Recently several authors ^ have proposed a generalized simulated 
annealing method to locate the global minimum based on the generalized 
statistical mechanics proposed by Tsallis [^, which can overcome the slow 
convergence of the traditional simulated annealing method 0] based on the 
standard Monte Carlo scheme. Subsequently, Andricioaei and Straub have 
pointed out that their generalized Monte Carlo (Tsallis Monte Carlo) |^] 
method may be used to overcome the quasi-ergodicity appearing in the 
canonical average at constant temperatures. 

In this short note we use a new Monte Carlo scheme ^, 10 1 based 
on Tsallis's generalized statistical mechanics |^] for the calculation of the 
thermal average of the potential energy of a simple double well potential at 
constant temperatures. We examine the performance of this new algorithm 
designed to overcome the quasi-ergodicity. 



2 A generalized Monte Carlo scheme 

In the generalized statistical mechanics proposed by Tsallis p , a crucial role 
is played by the generalized entropy Sq defined as 

Sq = k ^ (1) 

q — I 

where g is a real number which characterizes the statistical mechanics, and 
Pi is the probability of states i. This entropy Sq becomes the usual Gibbs- 
Shannon entropy 5i = —k^pilnpi when q ^ 1. By maximizing the gen- 
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eralized entropy with the constraints 



y^^Pi = 1 and ^pfcj = const 



where is the energy spectrum, the statistical weight for the generaUzed 
canonical ensemble characterized by the parameter q is given by the Tsallis 
weight 



pj = exp (-/?e,) /Z^^ (3) 



with 



= =: r hi 



1 - (1 - q)l3ei (4) 



where li is the generalized energy and (3 is the Lagrange multiplier, which 
plays the role of the generalized temperature in the Tsallis statistical me- 
chanics. We should note here that the statistical weight is not given by the 
probability distribution pi but rather by pj in order to calculate thermal 
averages, as can be seen in (^). The generalized canonical partition function 
Zq in @ is given by 

Zg = Y^ exp (-^Ej/g) (5) 

In the limit g — > 1, ^ and (3 ^ f3 = 1/kT where T is the usual 
temperature , then the Tsallis weight (^) becomes the Boltamann weight 
Pi = exp(— /3ei)/Zi of the usual canonical ensemble. 

In this generalized statistical mechanics, Andricioaei and Straub Q have 
noted that since the average of an observable O is defined by 

<0>,= Y,pIO,, (6) 

the detailed balance condition should be written as 

p1W,,=p]Wj, (7) 

where Wij is the element of the transition matrix. Based on this observation, 
they Q pointed out that the generalized Monte Carlo algorithm originally 
proposed by Penna Q in his generalized simulated annealing method cannot 
satisfy the detailed balance. They |^ then proposed a new generalized Monte 
Carlo scheme, where the acceptance probability p of the Monte Carlo move 
is given by 

p = min l,exp(— /3Ae) (8) 
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which has a similar form to the familiar Metropolis algorithm 
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j» = min [1, exp(— /3Ae)] (9) 

where Ae is the increase of the generahzed energy while Ae is that of 
the usual energy. Andricioaei and Straub |^ combined this Monte Carlo 
algorithm with the usual simulated annealing technique Q to optimize the 
conformation of tetrapeptides. They showed that their algorithm is more 
effective than the standard simulated annealing based on the usual molecular 
dynamics or Monte Carlo methods. 

Recently, Andricioaei and Straub ^ pointed out that their generalized 
Monte Carlo scheme Q for simulated annealing can also be used to calculate 
the standard canonical ensemble averages at constant temperatures 

< o >,= ^O'^-f-ff (10) 
using the general rule of the importance sampling Q 

</>=E/^ = E(f)ff^ (11) 

Then, equation ( [lO|) can be calculated using the generalized Monte Carlo 
scheme from 

J2 Oi exp (^/3ei - exp (^-fSli^ 

< O >i = J — J _ ^ (12) 

Eexp i^- yPei - I3eijj exp(-/3ei) 



<Oexp(- (/3e-/?e)) >g 
< exp(- (/?e-^)) >q 



(13) 



where the generalized average <>q is defined by (|§). In practice, the equilib- 
rium thermal average in the standard canonical ensemble can be calculated 
by conducting these "generalized" Monte Carlo steps and accumulating the 
samples by 

^^^^^EO.exp( (»,-^.)) ^^^^ 

Eexp i^- i^Pej - f^jjj 

where the sum j runs over to the generalized Tsallis Monte Carlo steps pro- 
duced by eq.(^) instead of the usual Metropolis algorithm (|9|). The equation 
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(|T^ is a generalized expression, and we may choose (3 = (3. Then, the 
equation (|^) reduces to 

<o>.= ^"'''-'^fe:^''> . (15) 

X;exp(-/3(ej -€j)) 

which has been used in g, |]. Fr om now on, the actual calculation will be 
done by choosing (3 = (3. 

We would like to point out here that a similar reweighting technique is 
used in the histogram method [13| and the multicanonical ensemble method 1 14] 



Andricioaei and Straub applied this Tsallis Monte Carlo algorithm to the 
two-dimensional Ising system |^], a classical one-dimensional potential and 
a 13-atom Lennard- Jones cluster Q , and showed that this new Monte Carlo 
algorithm is more effective than the standard Metropolis algorithm. They 
also proposed the generalized molecular dynamics based on Tsallis general- 
ized statistical mechanics 10 1 



3 Application to a classical double well potential 

In order to check the performance of this new Monte Carlo algorithm ^] 
in a classical system such as molecules and liquids, we look at the problem 
of a classical particle in a one-dimensional double well potential defined by 

y(x) = 5 (Sx*^ + 4(a - l)x^ - Gqx^) + 1 (16) 

where 

S = (17) 

2a + 1 ^ ' 

which was examined by Frantz et al. Q to demonstrate the quasi-ergodicity 
of the standard Metropolis algorithm. The potential given by ([l^ ) has the 
double minimum at x = 1 and x = —a; therefore, it represents a symmetrical 
double well when a = 1 and a single well when a = 0. This potential is 
more conveniently characterized by the parameter 7 Q 

^ F(0)-y(l) \2a + l) ^ ^ 

In figure 1, we show the asymmetric double well potential V{x) when 7 = 0.9 
as well as the equilibrium Tsallis weight oc exp(— /?y(x)), which is renor- 
malized, and Boltzmann weight pB = exp(— /?F(x))/Zi. It is obvious that 
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the statistical weight pT of the Tsalhs generahzed statistical mechanics has 
a greater chance of crossing the barrier and, therefore, a greater possibility 
of avoiding quasi-ergodicity. 

Alternatively this Tsallis's statistical mechanics can be regarded as a 
search-space smoothing method which deforms the rugged potential land- 
scape by a smoother one. In figure 2, we show the smoothed potential V{x) 
calculated from (^), which has a smoother landscape and a lower energy 
barrier than the original landscape V{x). 

We use this double well potential to look at the quasi-ergodicity of the 
standard as well as the generalized Monte Carlo scheme by examining the 
classical average potential energy 

^ ^ jV{x)ew{-(3V{x))dx ^^g^ 
/ exp (—(3V{x)) dx 

which can be calculated rigorously by numerical quadrature as a function of 
the temperature P and the asymmetry parameter 7. We employ the same 
Monte Carlo procedure of Frantz et al. and use the uniform sampling 
distribution from previous position x to new position x given by 

T(x \x) = — for la; — X | < — 

= otherwise (20) 

The step size A can be chosen to be large enough for the one-dimensional 
problem so that the quasi-ergodicity can be reduced However, we follow 
Frantz et al. and use the scaling A = 2.5/y^, which maintains an ap- 
proximately 50% acceptance rate at all temperatures (3, by considering the 
application to multidimensional systems. This 50% acceptance is commonly 
employed in the classical Monte Carlo community [15|. 



Since the error caused by the quasi-ergodicity depends upon the walk 
initialization, we start the Monte Carlo walk at the global minimum x = 1 
and at the metastable minimum x = —a. If the walk starts at the global 
minimum, then the average potential <V> will be low if the walk is quasi- 
ergodic, as the distribution associated with the higher energy well around 
a-minimum will be insufficiently sampled. While the walk is initialized at 
—a, the walk may be trapped at this metastable a-well and the average 
potential < V > will be too high. 

In figure 3, we show the average potential energy < V > calculated from 
(i) the standard Metropolis Monte Carlo and (ii~iv) the generalized Tsallis 
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Monte Carlo algorithm as a function of temperature (3. Similar curves can be 
found in figure 2 of reference Q; however, those authors plotted the results 
from the usual Metropolis algorithm initialized at the global minimum x = 1 
and at the random positions. In order to check the quasi-ergodicity of the 
walk, we instead show two curves initialized at the global minimum x = 1 
and the metastable a-minimum at x = —a for each scheme. The average 
potential < F > is obtained from 100 independently initialized walks, each 
consisting of 500 warm-up steps followed by 10^ steps with data accumu- 
lation. The average energy < V > calculated from the standard Metropo- 
lis algorithm shows large errors due to the quasi-ergodicity. In particular, 
the walk initialized at the metastable a-minimum greatly overestimates the 
potential energy. The walk starting from the global minimum also shows 
noticeable error and underestimates the average energy, as expected. 

In contrast to the results from the standard Metropolis algorithm, the 
results from the generalized Monte Carlo algorithm |9| show excellent 
agreement with the numerically exact result. The walks started from the 
global minimum x = 1 produce the average energy which is almost indis- 
tinguishable from the exact result up to the lowest temperature (high /?) 
examined. Even the walk started from the metastable a minimum produces 
the energies that are fairly accurate up to rather low temperature /? ~ 15 
for q = 1.5, which shows that this generalized algorithm does a much better 
job than the standard Monte Carlo scheme. The agreement with the exact 
result is further improved if we increase the value of parameter q from 1.5 
(see Figure 3). It should be noted that the probability distribution associ- 
ated with the parameter g = 2 is given by the Cauchy distribution and with 
g > 5/3 by the super-diffusive Levy distribution The theoretical upper 
boundary is g = 3. The quasi-ergodicity can be largely circumvented by 
the use of Tsallis's generalized statistical mechanics with q > 1, as shown in 
figure 3. 

In order to overcome the quasi-ergodicity, the Monte Carlo walk has to 
pass many times through the energy barrier. Figure 4 compares the time se- 
ries of the Monte Carlo walk generated by the standard Monte Carlo scheme 
{q = 1) and that of the generalized Monte Carlo scheme. Similar diagrams 
can be found in reference Q. Although the trajectory of the standard 
Monte Carlo scheme is mostly confined within the metastable a minimum, 
the number of barrier crossing of the Tsallis Monte Carlo scheme is signifi- 
cantly greater. This barrier crossing occurs more often when the parameter 
q is larger. The quasi-ergodicity can be removed by Tsallis's generalized 
Monte Carlo scheme by increasing the magnitude of the parameter q from 
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1. 

In figure 5, we show the average energy < V > as a function of the asym- 
metry parameter 7. Again a similar diagram can be found in reference [^. 
Here again we can clearly observe a large error for the initialization at the 
metastable a minimum when the standard Monte Carlo procedure is used. 
In this case, the superiority of Tsallis's generalized Monte Carlo scheme over 
the traditional one is also obvious. 

In table 1 we showed the effect of quasi-ergodicity in the convergence 
when 7 = 0.9 and (3 = 10 for both the standard Monte Carlo and the general- 
ized Monte Carlo algorithms with q = 1.5, each originating at the metastable 
a-minimum. We show the average energy <V> and the standard deviation 
(STD) obtained from 100 independently initialized walks. We found that 
the problem of quasi-ergodicity obtained with the standard Metropolis algo- 
rithm is still present, even if we increase the number of Monte Carol steps. 
On the other hand the result from the generalized Monte Carlo simulation 
is significantly improved by increasing the number of Monte Carlo steps, 
and ergodicity can be recovered. Furthermore, STD decreases as ~ 
when the Monte Carlo step N is increased for the generalized Monte Carlo 
simulation, while such a systemic decrease is not observed for the standard 
Monte Carlo simulation. 

This convergence can be clearly visualized by plotting the visiting prob- 
ability P{x) of the Monte Carlo walks, which should be given by the Boltz- 
mann weight for the standard Monte Carlo scheme and by the Tsallis weight 
@ for the generalized Monte Carlo scheme. In figure 6, we compare P{x) 
generated by the standard and generalized Monte Carlo walks with the 
Boltzmann and the Tsallis weights. The visiting probability P{x) converges 
rapidly to the Tsallis weight as the number of the Mone Carlo walks is in- 
creased for the generalized Monte Carlo scheme. On the other hand, P{x) 
recovers the Boltzmann weight only when the Monte Carlo steps are in- 
creased upto ~ 10^ in the standard Monte Carlo scheme. 

It is also interesting to discuss the convergence in the context of the " gen- 
eralized ergodic measure" of Thirumalai, Mountain, and Kirkpatrick [16, ^. 
Using the energy metric, the ergodic measure dY{n) is defined by 

dv{n) = {y%n)-V\n))\ (21) 

where V°'{n) and V^{n) are the move average of the potential energy V 
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along the trajectories a and b, which are defined from (IT5|) by 



y"(n) = ' ' / / ' (22) 

E".oexp(-/3(l^«-F;)) 

for the generahzed Tsalhs Monte Carlo scheme, and 

V-in) = ' ' ' ) il. (23) 

for the standard Monte Carlo scheme. For an ergodic system, Thirumalai 
et al. [16 1 suggested that the ergodic measure converges as dyiji) — > oo if 



n — > oo. They found that 



dy{n)^dv{^)^ (24) 



where the "diffusion" constant Dy depends on temperature. In figure 7, 
we show dy(0)/(iy(n) for 7 = 0.9 potential obtained from 100 independent 
pairs (a, b) of Monte Carlo walks as a function of Monte Carlo steps n. We 
chose the trajectory starting from the metastable minimum x = —a (=- 
0.9163) and that from the stable minimum x = 1, respectively, as the two 
independent trajectories a and b. It is clear from figure 7 that the normalized 
inverse of the ergodic measure dy{Q) / dy{n) grows linearly with the Monte 
Carlo steps n. Therefore, the ergodic measure dy{n) decreases according to 
(p^, and the diffusion constant Dy can be determined. 

Figure 8 shows the diffusion constant Dy as a function of temperature (3. 
The diffusion constant is a decreasing function of the inverse temperature 
/3; it is more difficult to recover the ergodicity when the temperature is 
lowered. However, at lower temperatures /? > 20 with f3/S.V > 1 where 
AV = V{—a) — V{1) = 0.1, the diffusion constant starts to increase again 
for the Tsallis Monte Carlo walks {q ^ 1). This is due to the fact that the 
equilibrium thermal distribution around the metastable minimum at x = —a 
becomes negligibly small at such low temperatures, and the ergodicity of the 
trajectory among the stable and the metastable minima is irrelevant once 
trajectories a and b both fall into the stable minimum around x = 1. 

It is apparent from figure 8 that the diffusion constant Dy obtained from 
the standard Monte Carlo algorithm obeys the usual activation form 

Dy ~ exp i-PEh) (25) 
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where the activation energy Ef, is estimated to be Ef, ~ 0.4, which is the same 
order of magnitude as the barrier height (~ 1) between two weUs, as shown 
in figure 1. On the other hand, the diffusion constant Dy obtained from 
the Tsalhs generalized Monte Carlo scheme is characterized by a power-law 
form 

Dv ~ (26) 

with the exponent 7]. This result is conceivable, based on the acceptance 
probability (^) using the Tsallis weight, because 

Dv ~ exp(— /3£'fe) 

~ [l-{l-q)PE,]-'^/^'^-'^ (27) 

and Dy ~ We note in figure 8 that rj ~ 1.6 for q = 1.5 and 

rj ~ 0.85 for q = 2.5. These values are close to q/{q — 1) = 3 for q = 
1.5 and q/{q — 1) = 1.66 if q = 2.5. Therefore, the diffusion constant 
decreases more slowly as a function of the inverse temperature /? as for 
the generalized Monte Carlo scheme, though it decreases exponentially faster 
for the standard Monte Carlo scheme. The exponent ry becomes smaller as 
q is increased and the diffusion is possible even at low temperatures. This 
result can be anticipated based on data presented in figures 3 and 6. This 
power-law temperature dependence of the diffusion constant Dy for g 7^ 1 
has been pointed out by Straub and Andricioaei |17|. 



4 Concluding remarks 

In conclusion, we found that the so-called "quasi-ergodicity" in a double well 
potential encountered in the standard Metropolis Monte Carlo algorithm is 
largely circumvented by the use of the generalized Monte Carlo algorithm 
of Andricioaei and Straub [^, ^ based on Tsallis's [Q] generalized statistical 
mechanics. Therefore this algorithm will be useful in the Monte Carlo study 
of any system, in particular, when the first order phase transition occurs. 
We believe that this algorithm will be useful, for example, in the study of 
the melting transition of clusters for which the J-walk algorithm M and the 



multiple histrogram method |18] have been utilized. It is also interesting 
to apply this new algorithm to the problems where the slow dynamics due 
to randomness or frustration is serious, for example, the spin-glass prob- 
lems. This generalized Monte Carlo algorithm based on Tsallis's statisti- 
cal mechanics use the reweighting technique similar to the multicanonical 
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method |14]. But this algorithm is much simpler than the multicanonical 



methods and appears to be more promising because it does not require an 
iterative determination of the weighting factor by preliminary Monte Carlo 
runs. 

Recently, a new formulation of the Tsallis statistics has been proposed |19] 



With a new choice of energy constraint, unfamiliar consequences in the pre- 
vious formulation have been erased. However, the Lagrange parameter used 
in the new formulation can be related to that in the previous formulation. 
Therefore, in performing Monte Carlo simulations with reweighting tech- 
nique, we may use the previous formulation of the Tsallis statistics. 
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Figure 1: A classical double well potential (16) (thin solid curve),for 7 = 0.9, 
the Tsallis weight oc exp{—pv{x)) with q = 1.5 (dotted curve), q = 2.5 
(chain curve) and the Boltzmann weight ps = exp{— PV (x)) / Zi (thick solid 
curve) both at temperature (3 = 10. The Tsallis weights are renormalized. 
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Figure 2: A smoothed potential landscape V{x) when q = 1.5 (dotted curve) 
and q = 2.5 (chain curve) for 7 = 0.9 and fi = 10, compared with the original 
potential landscape V{x) {q = 1, thin solid curve). The energy barrier is 
lowered by the smoothing defined by (3) of Tsallis's statistical mechanics. 
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Figure 3: The classical average potential energy < V > for the 7 = 0.9 
potential as a function of /3. The thin solid curves are initialized at the 
metastable a-minimum xq = —0.9613 with the standard Monte Carlo algo- 
rithm with the parameter q = 1 and the generalized Monte Carlo algorithm 
with q = 1.5, q = 2.0 and q = 2.5. The exact result is represented by the 
thick solid curve. The results initialized at the global minimum xq = 1 
and calculated from the Metropolis algorithm with q = 1 and that from the 
generalized Monte Carlo algorithm with q = 1.5 arc shown by the dotted 
curve. All results are obtained from 100 independently initialized Monte 
Carlo walks, each consisting of 500 warm-up steps followed by 10^ steps 
with data accumulation. Some representative error bars are also shown. 
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Figure 4: Time series of the Monte Carlo walks generated by (a) the 
standard Monte Cairo (q = 1), (b) th$()rsallis Monte Carlo with q = 1.5, and 
(c) q = 2.5 when /3 = 10 and 7 = 0.9. All walks start from the metastable a- 
minimum, and 500 warm up steps are also included. The number of barrier 
crossing is greater for larger q. 
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Figure 5: The average potential energy < F > for the temperature /3 = 10 
as a function of the asymmetry parameter 7 of the double well potential. 
The solid and dotted curves are used the same as in figure 3. We only show 
the result and some representative error bars for the cases q = 1 and q = 1.5 
of figure 3. 
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Figure 6: The Visiting probabilities P{x) obtained from (a) the standard 
Monte Carlo walk compared with the Boltzmann weight, and (b) the gener- 
alized Monte Carlo walks compared with the Tsallis weight (c.f. figure 1). 
Probability distribution P{x) generated from the generalized Monte Carlo 
scheme can recovers the equilibrium Tsallis weight within short Monte Carlo 
steps (a), while it cannot recover the equilibrium Boltzmann weight until 
rather long (~ 10^) Monte Carlo steps are executed (b). 
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MC steps 

Figure 7: The normalized inverse of the generalized crgodic measure 
dy {0) / dv {n) as a function of the number of Monte Carlo steps n calcu- 
lated for 7 = 0.9 potential at /? = 10 using the standard Monte Carlo scheme 
{q = 1) and the generalized Monte Carlo scheme {q = 1.5, 2.5). dv{0)/dv{n) 
grows linearly with the Monte Carlo steps n and the diffusion is faster for 
larger q. 
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Figure 8: The diffusion constant Dy as a function of temperature (3. (a) A 
straight hne of the semi-log plot of Dy vs (3 for the standard Monte Carlo 
scheme indicates activation-type diffusion, while (b) straight lines of the log- 
log plot for the Tsallis Monte Carlo scheme in the high temperature region 
indicate slow power-law decrease of the diffusion constant as the temperature 
is lowered (larger j3). 
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Table 1: The convergence of the Metropohs and the generahzed Tsalhs 
Monte Carlo algorithms at /? = 10. The average energy <V>, the standard 
deviation (STD) obtained from 100 independently initialized walks, each 
originating at the metastable a-minimum. 





Metropolis Monte Carlo 


Tsallis Monte Carlo {q = 1.5) 


steps 


<V> (STD) 


<V> (STD) 


10^ 


0.1500 (0.0177) 


0.1065 (0.0518) 


10^ 


0.1520 (0.0049) 


0.0937 (0.0394) 


10^ 


0.1475 (0.0192) 


0.0798 (0.0162) 


10^ 


0.1052 (0.0320) 


0.0798 (0.0049) 


10^ 


0.0805 (0.0124) 


0.0800 (0.0017) 


Exact 


0.0799 


0.0799 
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